function Phi = basis_function_loc(V, x, y, basis_type, derivative_type)
Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
J_det = det(Jacobi);
x_hat = (Jacobi(2,2)*(x-V(1,1)) - Jacobi(1,2)*(y-V(1,2))) / J_det;
y_hat = (-Jacobi(2,1)*(x-V(1,1)) + Jacobi(1,1)*(y-V(1,2))) / J_det;
if strcmp(derivative_type, "x")
    Phi = basis_function_ref(x_hat, y_hat, basis_type, "x");
else
    dNdx = basis_function_ref(x_hat, y_hat, basis_type, "dx");
    dNdy = basis_function_ref(x_hat, y_hat, basis_type, "dy");
    if strcmp(derivative_type, "dx")
        Phi = (Jacobi(2,2) * dNdx - Jacobi(2,1) * dNdy) / J_det;
    elseif strcmp(derivative_type, "dy")
        Phi = (-Jacobi(1,2) * dNdx + Jacobi(1,1) * dNdy) / J_det;
    else
        error("Invalid derivative type.");
    end
end
end

%% basis_function_ref
function N = basis_function_ref(x, y, basis_type, derivative_type)
if strcmp(basis_type, "P1")
    N = zeros(size(x,1), 3);
elseif strcmp(basis_type, "P1b")
    N = zeros(size(x,1), 4);
elseif strcmp(basis_type, "P2")
    N = zeros(size(x,1), 6);
else
    error("Invalid basis type.");
end
if strcmp(basis_type, "P1")
    if strcmp(derivative_type, "x")
        N(:,1) = 1 - x - y;
        N(:,2) = x;
        N(:,3) = y;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
    end
elseif strcmp(basis_type, "P1b")
    if strcmp(derivative_type, "x")
        N(:,1) = 1 - x - y;
        N(:,2) = x;
        N(:,3) = y;
        N(:,4) = 27*x.*y.*N(:,1);
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
        N(:,4) = -27*y.*(2*x + y - 1);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
        N(:,4) = -27*x.*(x + 2*y - 1);
    end
elseif strcmp(basis_type, "P2")
    lambda = 1 - x - y;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda - 1);
        N(:,4) = 4*x.*lambda;
        N(:,2) = x.*(2*x - 1);
        N(:,5) = 4*x.*y;
        N(:,3) = y.*(2*y - 1);
        N(:,6) = 4*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 1 - 4*lambda;
        N(:,4) = 4*(lambda - x);
        N(:,2) = 4*x - 1;
        N(:,5) = 4*y;
        N(:,3) = 0;
        N(:,6) = -4*y;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 1 - 4*lambda;
        N(:,4) = -4*x;
        N(:,2) = 0;
        N(:,5) = 4*x;
        N(:,3) = 4*y - 1;
        N(:,6) = 4*(lambda - y);
    end
end
end